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Abstract — Wavelets have been used extensively for several 
years now in astronomy for many purposes, ranging from data 
filtering and deconvolution, to star and galaxy detection or 
cosmic ray removal. More recent sparse representations such 
ridgelets or curvelets have also been proposed for the detection of 
anisotropic features such cosmic strings in the cosmic microwave 
background. We review in this paper a range of methods based on 
sparsity that have been proposed for astronomical data analysis. 
We also discuss what is the impact of Compressed Sensing, 
the new sampling theory, in astronomy for collecting the data, 
transferring them to the earth or reconstructing an image from 
incomplete measurements. 

Index Terms — Astronomical data analysis. Wavelet, Curvelet, 
restoration, compressed sensing 



I. Introduction 

The wavelet transform (WT) has been extensively used 
in astronomical data analysis during the last ten years. A 
quick search with ADS (NASA Astrophysics Data System, 
adswww.harvard.edu) shows that around 1000 papers contain 
the keyword "wavelet" in their abstract, and this holds for 
all astrophysical domains, from study of the sun through to 
CMB (Cosmic Microwave Background) analysis ||29[] . This 
broad success of the wavelet transform is due to the fact that 
astronomical data generally gives rise to complex hierarchi- 
cal structures, often described as fractals. Using multiscale 
approaches such as the wavelet transform, an image can 
be decomposed into components at different scales, and the 
wavelet transform is therefore well-adapted to the study of 
astronomical data. Furthermore, since noise in the physical 
sciences is often not Gaussian, modeling in wavelet space of 
many kind of noise - Poisson noise, combination of Gaussian 
and Poisson noise components, non-stationary noise, and so 
on - has been a key motivation for the use of wavelets in 
astrophysics. 

If wavelets represent well isotropic features, they are far 
from optimal for analyzing anisotropic objects. This has 
motived other constructions such the curvelet transform 
More generally, the best data decomposition is the one which 
leads to the sparsest representation, i.e. few coefficients have a 
large magnitude, while most of them are close to zero. Hence, 
for specific astronomical data set containing edges (planetary 
images, cosmic strings, etc), curvelets should be preferred to 
wavelets. 
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Fig. 1. Galaxy NGC 2997. 



In this paper, we review a range of astronomical data 
analysis methods based on sparse representations. We first in- 
troduce the Undecimated Isotropic Wavelet Transform (UIWT) 
which is the most popular WT algorithm in astronomy. We 
show how the signal of interest can be detected in wavelet 
space using a noise modeling, allowing us to build the so- 
called multiresolution support. Then we present inlllllhow this 
multiresolution support can used for restoration applications. 
In section IIVI another representation, the curvelet transform, 
is introduced, which is well adapted to anisotropic structure 
analysis. Combined together, the wavelet and the curvelet 
transforms are very powerful to detect and discriminate very 
faint features. We give an example of application for cosmic 
string detection. Section |V] describes the compressed sensing 
theory which is strongly related to sparsity, and presents its 
impacts in astronomy, especially for spatial data compression. 

II. The Isotropic Undecimated Wavelet Transform 



The Isotropic undecimated wavelet transform (lUWT) |i25| 
decomposes an n x n image cq into a coefficient set W = 
{wi, . . . , wj, cj}, as a superposition of the form 

,7 

co[k,l] = cj[k,l] +y^Wj[fc,Z], 

where cj is a coarse or smooth version of the original image 
Co and Wj represents 'the details of cq' at scale 2^^ (see 
Starck et al.[30, 28] for more information). Thus, the algorithm 
outputs J + 1 sub-band arrays of size n x n. (The present 
indexing is such that j = 1 corresponds to the finest scale 
(high frequencies)). 

Hence, we have a multi-scale pixel representation, i.e. each 
pixel of the input image is associated with a set of pixels 
of the multi-scale transform. This wavelet transform is very 
well adapted to the detection of isotropic features, and this 
explains its success for astronomical image processing, where 
the data contain mostly isotropic or quasi-isotropic objects, 
such as stars, galaxies or galaxy clusters. 

The decomposition is achieved using the filter bank 
(^2D, g2D = S~h2D, h2D = S, g2D = S) whcrc h2D is the tensor 
product of two ID filters hn). The passage from one resolution 
to the next one is obtained using the "a trous" algorithm I30ll 



Cj+i[k,l] = hiD[m]hiD[n]cj[k + 2^m, I + 2^n] 

m n 

Wj+i[fc, /] = Cj[kJ] - Cj+i[kJ] , 



(1) 
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Fig. 2. Wavelet transform of NGC 2997 by the lUWT. The co-addition of these six images reproduces exactly the original image. 



where huj is typically a symmetric low-pass filter such as the 
£?3 -Spline filter 

Fig. |2] shows lUWT of the galaxy NGC 2997 displayed in 
Fig. [1] . Five wavelet scales are shown and the final smoothed 
plane (lower right). The original image is given exactly by the 
sum of these six images. 

A. Example: Dynamic range compression using the lUWt 



Fig. 3. Top left - Hale-Bopp Comet image. Bottom left - histogram 
equalization results. Bottom right - wavelet-log representations. 

Since some features in an image may be hard to detect 
by the human eye due to low contrast, we often process the 
image before visualization. Histogram equalization is certainly 
one the most well-known methods for contrast enhancement. 
Images with a high dynamic range are also difficult to analyze. 
For example, astronomers generally visualize their images 
using a logarithmic look-up-table conversion. 

Wavelets can be used to compress the dynamic range at all 
scales, and therefore allow us to clearly see some very faint 
features. For instance, the wavelet-log representation consists 
of replacing Wj[k, I] by sgn{wj[k, I]) log{\wj[k, l]\), leading to 
the alternative image 

,7 

Ik,i = log(cj,A;,;) + ^sgn(wj[A:,;])log(| Wj[k,l] \ +e) (2) 

where e is a small number (for example e = 10^'^). Fig. [3] 
shows a Hale-Bopp Comet image (logarithmic representation) 
(top left), its histogram equalization (middle row), and its 
wavelet-log representation (bottom). Jets clearly appear in the 
last representation of the Hale-Bopp Comet image. 

B. Signal detection in the wavelet space 

Observed data Y in the physical sciences are generally 
corrupted by noise, which is often additive and which follows 
in many cases a Gaussian distribution, a Poisson distribution, 
or a combination of both. It is important to detect the wavelet 
coefficients which are "significant", i.e. the wavelet coeffi- 
cients which have an absolute value too large to be due to 
noise. We defined the multiresolution 1\I of an image Y by: 

1 if Wj[k^V\ is significant 
if Wj[k^V\ is not significant 

where Wj [fc, I] is the wavelet coefficient of Y at scale j and at 
position (fc,Z). We need now to determine when a wavelet 
coefficient is significant. For Gaussian noise, it is easy to 
derive an estimation of the noise standard deviation cr-j at scale 



Mj[k,l] 



(3) 



i from the noise standard deviation, which can be evaluated 
with good accuracy in an automated way 02711 . To detect 
the significant wavelet coefficients, it suffices to compare the 



wavelet coefficients Wjlk^l] to a threshold level tj. tj is 
generally taken equal to Kaj, and K is chosen between 3 
and 5. The value of 3 corresponds to a probability of false 
detection of 0.27%. If wj [k, I] is small, then it is not significant 
and could be due to noise. If Wj[k, I] is large, it is significant: 



if I Wj[k,l] I > tj then Wj[k,l] is significant 
if I Wj[k,l] I < tj then Wj[fc,Z] is not significant 



(4) 



When the noise is not Gaussian, other strategies may be 
used: 

• Poisson noise: if the noise in the data Y is Poisson, 
the ti-ansformation Jl] A{Y) = '^\J I + | acts as if the 
data arose from a Gaussian white noise model, with 
(7 = 1, under the assumption that the mean value of / 
is sufficiently large. However, this transform has some 
limits and it has been shown that it cannot be applied for 
data with less than 20 photons per pixel. So for X-ray 
or gamma ray data, other solutions have to be chosen, 
which manage the case of a reduced number of events or 
photons under assumptions of Poisson statistics 

• Gaussian + Poisson noise: the generalization of variance 
stabilization |18] is: 



aY[k,l] 



ag 



where a is the gain of the detector, and g and a are the 
mean and the standard deviation of the read-out noise. 
Poisson noise with few events using tlie MS-VST 

For images with very few photons, one solution consists 
in using the Multi-Scale Variance Stabilization Trans- 
form (MS VST) 132]. The MS VST combines both the 
Anscombe transform and the lUWT in order to produce 
stabilized wavelet coefficients, i.e. coefficients corrupted 
by a Gaussian noise with a standard deviation equal to 1 . 
In this framework, wavelet cofficients are now calculated 
by: 

lUWT r c, = Y.mY.r^hiD[m\h^D[n] 

+ I cj_i[fc + 2^-^771,; + 2^-171] (5) 

MS-VST [ W.J = Aj-iicj-i) -Ajicj) 

where Aj is the VST operator at scale j defined by: 



A,ic,)^b^^^J\c,+eU)\ 



(6) 



where the variance stabilization constants b^-'-' and e^^^ 
only depends on the filter /iijj and the scale level j. They 
can all be pre-computed once for any given h [32]. The 
multiresolution support is computed from the MSVST 
coefficients, considering a Gaussian noise with a standard 
deviation equal to 1. This stabiUzation procedure is also 
invertible as we have: 



Co = ^0 



(7) 
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For other kind of noise (correlated noise, non stationary 
noise, etc), other solutions have been proposed to derive the 
multiresolution support |29]. In next section, we show how 
the multiresolution support can be used for denoising and 
deconvolution. 

III. Restoration using the Wavelet Transform 



Example 

A simulated Hubble Space Telescope image of a distant 
cluster of galaxies is shown in Fig. ID middle. The simulated 
data are shown in Fig. HI left. Wavelet deconvolution solution 
is shown Fig. |4] right. The method is stable for any kind of 
point spread function, and any kind of noise modeling can be 
considered. 



A. Denoising 

The most used filtering method is the hard thresholding, 
which consists of setting to all wavelet coefficients of Y 
which have an absolute value lower than a threshold i„ 



Wj [k, I] 







,-[fc, I] if I Wj[k, I] |> tj 
otherwise 



(8) 



More generally, for a given sparse representation (wavelet, 
curvelet, etc) with its associate fast transform Tu, and fast 
reconstruction TZ^, we can derive a hard thresholding de- 
noising solution X from the data Y , by first estimating the 
multiresolution support M using a given noise model, and 
then calculating: 

X = n^M%,Y. (9) 

We transform the data, multiply the coefficients by the support 
and reconstruct the solution. 

The solution can however be improved considering the 
following optimization problem minx || M{%uY — %uX) jjj 
where M is the multiresolution support of Y. A solution can 
be obtained using the Landweber iterative scheme |22, 30ll : 

n^,M[Tu,Y-%,X"] (10) 



If the solution is known to be positive, the positivity constraint 
can be introduced using the following equation: 



- p+ (X" + 7^,„M [%,Y - r^x"] 



(11) 



where P+ is the projection on the cone of non-negative images. 

This algorithm allows us to constraint the residual to have 
a zero value inside the multiresolution support [,30ll . For 
astronomical image filtering, iterating improves significantly 
the results, especially for the photometry (i.e. the integrated 
number of photons in a given object). 



B. Deconvolution 

In a deconvolution problem, Y — HX+N, when the sensor 
is linear, H is the block Toeplitz matrix. Similarly to the 
denoising problem, the solution can be obtained minimizing 
minx II A/7^(y — HX) II2 under a positivity constraint, 
leading to the Landweber iterative scheme 11221 BOll : 



^n+i ^ i^j^u ^ H^Tl^MT^ [Y - HX"]) 



(12) 



Only coefficients that belong to the multiresolution support are 
kept, while the others are set to zero 12211 . At each iteration, the 
multiresolution support AI can be updated by selecting new 
coefficients in the wavelet transform of the residual which have 
an absolute value larger than a given threshold. 



C. Inpainting 

Missing data are a standard problem in astronomy. They can 
be due to bad pixels, or image area we consider as problematic 
due to calibration or observational problems. These masked 
area lead to many difficulties for post-processing, especially 
to estimate statistical information such the power spectrum or 
the bispectrum. The inpainting technique consists in filling the 
gaps. The classical image inpainting problem can be defined 
as follows. Let X be the ideal complete image, Y the observed 
incomplete image and L the binary mask (i.e. L[k,l] = 1 if 
we have information at pixel (fc, I), L[k, I] = otherwise). In 
short, we have: Y = LX. Inpainting consists in recovering X 
knowing Y and L. 

Noting llzllo the Iq pseudo-norm, i.e. the number of non- 
zero entries in z and ||z|| the classical I2 norm (i.e. ||z|p = 
Sfcl-^fc)^)' ^^^^ want to minimize: 

$'^X| 



mm I 
X 



subject to II y - LX \\i,^< cr, (13) 



where a stands for the noise standard deviation in the noisy 
case. It has also been shown that if X is sparse enough, 
the Iq pseudo-norm can also be replaced by the convex li 
norm (i.e. ||z||i = X]fc l^fcD lE3l- The solution of such 
an optimization task can be obtained through an iterative 
thresholding algorithm called MCA Q : 



X 



ri+l 



Y -LX") 



(14) 



where the nonlinear operator A$.a(-Z^) consists in: 

> decomposing the signal Z on the dictionary $ to derive 

the coefficients a = 
• threshold the coefficients: a = p{a, A), where the thresh- 
olding operator p can either be a hard thresholding (i.e. 
p{ai,X) — ai if \ai\ > A and otherwise) or a soft 
thresholding (i.e. p{ai,X) — sign(Q!i)max(0, \ai\ — A)). 
The hard thresholding corresponds to the Iq optimization 
problem while the soft-threshold solves that for li. 
m reconstruct Z from the thresholds coefficients a. 

The threshold parameter A„ decreases with the iteration num- 
ber and it plays a role similar to the cooling parameter of 
the simulated annealing techniques, i.e. it allows the solution 
to escape from local minima. More details relative to this 
optimization problem can be found in \12, [l^. For many 
dictionaries such as wavelets or Fourier, fast operators exist 
to decompose the signal so that the iteration of eq. [14] is very 
fast. It requires only to perform at each iteration a forward 
transform, a thresholding of the coefficients and an inverse 
transform. 

Example: The experiment was conducted on a simulated 
weak lensing mass map masked by a typical mask patterns 
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Fig. 4. Simulated Hubble Space Telescope image of a distant cluster of galaxies. Left: original, unaberrated and noise-free, middle: input, aberrated, noise 
added. Right, wavelet restoration wavelet. 



Fig. 5. Left panel, simulated weak lensing mass map, middle panel, simulated mass map with the mask pattern of CFHTLS data, right panels, inpainted 
mass map. The region shown is 1° x 1°. 



Fig. 6. A few first generation curvelets. 

(see Fig. |5]l. The left panel shows the simulated mass map 
and the middle panel show the masked map. The result of the 
inpainting method is shown in the right panel. We note that 
the gaps are undistinguishable by eye. More interesting, it has 
been shown that, using the inpainted map, we can reach an 
accuracy of about 1% for the power spectrum and 3% for the 
bispectrum 1.19,1 . 



IV. From Wavelet to Curvelet 

The 2D curvelet transform |9] was developed in an attempt 
to overcome some limitations inherent in former multiscale 
methods e.g. the 2D wavelet, when handling smooth images 
with edges i.e. singularities along smooth curves. Basically, 
the curvelet dictionary is a multiscale pyramid of localized 
directional functions with anisotropic support obeying a spe- 
cific parabolic scaling such that at scale 2^^, its length is 
2~^^^ and its width is 2~^. This is motivated by the parabolic 
scaling property of smooth curves. Other properties of the 
curvelet transform as well as decisive optimality results in 
approximation theory are reported in |8]. Notably, curvelets 
provide optimally sparse representations of manifolds which 
are smooth away from edge singularities along smooth curves. 
Several digital curvelet transforms lEsl ^ have been pro- 
posed which attempt to preserve the essential properties of 
the continuous curvelet transform and several papers report 
on their successful application in astrophysical experiments 

Fig. |6] shows a few curvelets at different scales, orientations 
and locations. 



Application to the detection of cosmic strings 

Some applications require the use of sophisticated statistical 
tools in order to detect a very faint signals, embedded in noise. 
An interesting case is the detection of non-Gaussian signatures 
in Cosmic Microwave Background (CMB). which is of great 
interest for cosmologists. Indeed, the non-Gaussian signatures 
in the CMB can be related to very fundamental questions 
such as the global topology of the universe [20], superstring 
theory, topological defects such as cosmic strings |6], and 
multi-field inflation 101 ■ The non-Gaussian signatures can, 
however, have a different but still cosmological origin. They 
can be associated with the Sunyaev-Zel'dovich (SZ) effect 13111 
(inverse Compton effect) of the hot and ionized intra-cluster 
gas of galaxy clusters yj], with the gravitational lensing by 
large scale structures, or with the reionization of the universe 



111]. They may also be simply due to foreground emission, or 
to non-Gaussian instrumental noise and systematics . 

All these sources of non-Gaussian signatures might have 
different origins and thus different statistical and morpholog- 
ical characteristics. It is therefore not surprising that a large 
number of studies have recently been devoted to the subject 
of the detection of non-Gaussian signatures. In H, 21], it was 
shown that the wavelet transform was a very powerful tool to 
detect the non-Gaussian signatures. Indeed, the excess kurtosis 
(4th moment) of the wavelet coefficients outperformed all the 
other methods (when the signal is characterized by a non-zero 
4th moment). 

Finally, a major issue of the non-Gaussian studies in CMB 
remains our ability to disentangle all the sources of non- 
Gaussianity from one another. It has been shown it was 
possible to separate the non-Gaussian signatures associated 
with topological defects (cosmic strings) from those due to the 
Doppler effect of moving clusters of galaxies (i.e. the kinetic 
Sunyaev-Zel'dovich effect), both dominated by a Gaussian 
CMB field, by combining the excess kurtosis derived from 
both the wavelet and the curvelet transforms 112 111 . 

The wavelet transform is suited to spherical-like sources 
of non-Gaussianity, and a curvelet transform is suited to 
structures representing sharp and elongated structures such as 
cosmic strings. The combination of these transforms highlights 
the presence of the cosmic strings in a mixture CMB-i-SZ-i-CS. 
Such a combination gives information about the nature of the 
non-Gaussian signals. The sensitivity of each transform to a 
part icular shape makes it a very strong discriminating tool 



Fig. 7. Top, primary Cosmic Microwave Background anisotropics (left) 
and kinetic Sunyaev-Zel'dovich fluctuations (right). Bottom, cosmic string 
simulated map (left) and simulated observation containing the previous three 
components (right). The wavelet function is overplotted on the Sunyaev- 
Zel'dovich map and the curvelet function is overplotted on the cosmic string 
map. 



In order to illustrate this, we show in Fig. [T] a set of 
simulated maps. Primary CMB, kinetic SZ and cosmic string 
maps are shown respectively in Fig. |7] top left, top right 
and bottom left. The "simulated observed map", containing 
the three previous components, is displayed in Fig. [T] bottom 
right. The primary CMB anisotropics dominate all the signals 
except at very high multipoles (very small angular scales). 
The wavelet function is overplotted on the kinetic Sunyaev- 
Zel'dovich map and the curvelet function is overplotted on 
cosmic string map. 
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V. Compressed Sensing 

A. Compressed Sensing in a nutshell 

Compressed sensing (CS) ifiol [Till is a new sam- 
pling/compression theory based on the revelation that one can 
exploit sparsity or compressibility when acquiring signals of 
general interest, and that one can design nonadaptive sampling 
techniques that condense the information in a compressible 
signal into a small amount of data. The gist of Compressed 
Sensing (CS) relies on two fundamental properties : 

1) Compressibility of the data : The signal X is said to 
be compressible if it exists a dictionary $ where the 
coefficients a = ^^X, obtained after decomposing X 
on <i>, are sparsely distributed. 

2) Acquiring incoherent measurements : In the Com- 
pressed Sensing framework, the signal X is not acquired 
directly; one then acquires a signal X by collecting data 
of the form Y — AX + rj : A is an m x n (with 
m < n) "sampling" or measurement matrix, and rj is 
a noise term. Assuming X to be sparse, the incoherence 
of A and $ (e.g. the Fourier basis and the Dirac basis) 
entails that the information carried by X is diluted in all 
the measurements Y. Combining the incoherence of A 
and $ with the sparsity of X in $ makes the decoding 
problem tractable. 

In the following, we choose the measurement matrix A to 
be a submatrix of an orthogonal matrix : the resulting 
measurement matrix is denoted Qa and obtained by picking 
a set of columns of Q indexed by A; 9a is obtained by 
subsampling the transformed signal QX. In practice, when Q 
admits a fast implicit transform (i.e. discrete Fourier transform, 
Hadamard transform, noiselet transform), the compression 
step is very fast and made reliable for on-board satellite 
implementation. 

A standard approach in CS attempts to reconstruct X by 
solving 

min \\a\\i^ s. t. ||y - QA^a\\e^ < e (15) 

a 

where is an estimated upper bound on the noise power. 

B. Compressed sensing for the Herschel data 

The Herschel/PACS mission of the European Space Agency 
(ESA) is facing with a strenuous compression dilemma : it 
needs a compression rate equal Xo p — 1/N with = 6. A 
first approach has been proposed which consists in averaging 
N = 6 consecutive images of a raster scan and transmitting 
the final average image. Nevertheless, doing so with high 
speed raster scanning leads to a dramatic loss in resolution. 
In jsl, we emphasized on the redundancy of raster scan 
data : 2 consecutive images are almost the same images up 
to a small shift 6. Then, jointly compressing/decompressing 
consecutive images of the same raster scan has been put 
forward to alleviate the Herchel/PACS compression dilemma. 
The problem then consists in recovering a single image X 
from N compressed and shifted noisy versions of X : 

yi e {!,■■■ ,N}; X,^TsAX) + m (16) 



where Ts- is an operator that shifts the original image X with 
a shift Si . The term 7/j models instrumental noise or model im- 
perfections. According to the compressed sensing framework, 
each signal is projected onto the subspace ranged by G. Each 
compressed observation is then obtained as follows : 

Vze {I,--- ,iV}; Y, = e,A,X, (17) 

where the sets {Aj} are such that the union of all the 
measurement matrices [8ai, • • • jQaJ span M". In practice, 
the subsets A, are disjoint and have a cardinality m = ln/N\. 
When there is no shift between consecutive images, these 
conditions guarantee that the signal X can be reconstructed 
univocally from {Yi}i=i^... at, up to noise. The decoding step 
amounts to seeking the signal x as follows : 

N 

min \\a\\i^ s. t. V \\Y, - eA.$a||fe < (18) 

i=l 

The solution of this optimization problem can be found via an 
iterative thresholding algorithm (see [5J) : 

N 

= A$^A„ (X" + tieJ2 - (19) 

1=1 

where the nonUnear operator A$ a(^) is defined in Equa- 
tion[T4]and the step-size /isa < "^/J^i \\'^Ai'^Ai\\2- Similarly 
to the MCA algorithm, the threshold A„ decreases with the 
iteration number towards the final value : A/; a typical value 
is A/ = 2 — 3(T. This algorithm has been shown to be very 
efficient for solving the problem in Equation [15] in 15|] . 



Fig. 8. Top left : Original image. Top right : First input noisy map (out of 6). 
The PACS data already contains approximately Gaussian noise. Bottom left : 
Mean of the 6 input images. Bottom right : Reconstruction from noiselet-based 
CS projections. The iterative algorithm has been used with 100 iterations. 



a) Illustration: We compare two approaches to solve 
the Herschel/PACS compression problem : i) transmitting the 
average of 6 consecutive images (M06), ii) compressing 6 
consecutive images of a raster scan and decompressing using 
Compressed Sensing. Real Herschel/PACS data are complex : 
the original datum X is contaminated with a slowly varying 
"flat field" component Cf. In a short sequence of 6 consecutive 
images, the flat field component is almost fixed. In this context, 
the data {xi}i=Q,... i can then be modeled as follows : 

X,^TsAX) + n, + Cf (20) 

If Cf is known (which will be the case in the forthcoming 
experiments), Tg. (X^"^) is replaced by Tg^ (AT^"^) + c/ in 
Equation [191 The data have been designed by adding realistic 
pointwise sources to real calibration measurements performed 
in mid-2007. In the following experiment, the sparsifying 
dictionary is $ is an undecimated wavelet tight frame and the 
measurement matrices are submatrices of the noiselet basis 
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The top-left picture of Figure [8] features the original signal X. 
In the top-right panel of Figure [8] The "flat field" component 
overwhelms the useful part of the data so that x has at best 
a level that is 30 times lower than the "flat field" component. 
The M06 solution (resp. the CS-based solution) is shown on 
the left (resp. right) and at the bottom of Figure [8] We showed 
in Jsll that Compressed Sensing provides a resolution enhance- 
ment that can reach 30% of the FWHM of the instrument's 
PSF for a wide range of signal intensities (i.e. flux of X). 
This experiment illustrates the reliability of the CS-based 
compression to deal with real-world data compression. The ef- 
ficiency of Compressed Sensing applied to the Herschel/PACS 
data compression relies also on the redundancy of the data : 
consecutive images of a raster scan are fairly shifted versions 
of a reference image. The good perfomances of CS is obtained 
by merging the information of consecutive images. The same 
data fusion scheme could be used to reconstruct with high 
accuracy wide sky areas from full raster scans. 

VI. Conclusion 

By establishing a direct link between sampling and sparsity, 
compressed sensing had a huge impact in many scientific 
fields, especially in astronomy. We have seen that CS could 
offer an elegant solution to the Herschel data transfer problem. 
By emphasing so rigorously the importance of sparsity, com- 
pressed sensing has also shed light on all work related to sparse 
data representation (such as the wavelet transform, curvelet 
transform, etc.). Indeed, a signal is generally not sparse in 
direct space (i.e. pixel space), but it can be very sparse after 
being decomposed on a specific set of functions. For inverse 
problems, compressed sensing gives a strong theoretical sup- 
port for methods which seek a sparse solution, since such a 
solution may be (under appropriate conditions) the exact one. 
Similar results are hardly accessible with other regularization 
methods. This explain why wavelets and curvelets are so 
successful for astronomical image denoising, deconvolution 
and inpainting. 
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